home *** CD-ROM | disk | FTP | other *** search
/ Precision Software Appli…tions Silver Collection 4 / Precision Software Applications Silver Collection Volume 4 (1993).iso / stats / chadyn.exe / YDBLRTR.C < prev    next >
C/C++ Source or Header  |  1988-11-28  |  5KB  |  163 lines
  1. /********************* YDblRtr.C for the double rotor ************************/
  2. /********************* (C) 1986,7,8 by JAMES A. YORKE ************************/
  3.  
  4.  
  5. #include "yinclud.h"
  6.  
  7. double  exp();
  8.  
  9. static double   u11,
  10.                 u12,
  11.                 u21,
  12.                 u22,
  13.                 s1,
  14.                 s2;
  15.  
  16.  
  17.  
  18. /* for double rotor: */
  19.  
  20. static double   Veloc11,
  21.                 Veloc12,
  22.                 Veloc21,
  23.                 Veloc22;
  24. static double   Matrix11,
  25.                 Matrix12,
  26.                 Matrix21,
  27.                 Matrix22;
  28. static double   cc1,
  29.                 cc2;
  30.  
  31.  
  32.  
  33. start33 () {            /* calculate parameters for the double rotor
  34.                    map */
  35.     double  T,
  36.             mass1,
  37.             mass2,
  38.             length1,
  39.             length2,
  40.             damping1,
  41.             damping2;
  42.  /* above parameters mostly to be set from menu */
  43.     double  inertia1,
  44.             inertia2;
  45.     double  discrim,
  46.             expon1,
  47.             expon2,
  48.             unorm;
  49.     mass1 = C1;        /* mass1 */
  50.     mass2 = C2;        /* mass2 */
  51.     length2 = C4;        /* length2 -- length1 is set below */
  52.     damping1 = C5;        /* damping1 & damping2 are damping  */
  53.     damping2 = C6;
  54.     T = C7;
  55.  
  56.  /* Let the length of the first rod and the inertias be set to values that
  57.     allow for a symmetric set of ODEs for the system */
  58.     length1 = length2 * sqrt(mass2 / (mass1 + mass2));
  59.     inertia1 = (mass1 + mass2) * length1 * length1;
  60.     inertia2 = mass2 * length2 * length2;
  61.  
  62.  /* compute eigenvalues for revised map */
  63.     discrim = 0.5 * sqrt(damping1 * damping1 + 4.0 * damping2 * damping2);
  64.     s1 = 0.5 * damping1 + damping2 + discrim;
  65.     s2 = 0.5 * damping1 + damping2 - discrim;
  66.     u11 = damping2 - s1;
  67.     u12 = damping2;
  68.     unorm = sqrt(u11 * u11 + u12 * u12);
  69.     u11 = u11 / unorm;
  70.     u12 = u12 / unorm;
  71.     u21 = damping2 - s2;
  72.     u22 = damping2;
  73.     unorm = sqrt(u21 * u21 + u22 * u22);
  74.     u21 = u21 / unorm;
  75.     u22 = u22 / unorm;
  76.  
  77.  /* Compute matrix M(T) for theta dot */
  78.     expon1 = exp(-s1 * T);
  79.     expon2 = exp(-s2 * T);
  80.     Veloc11 = u11 * u11 * expon1 + u21 * u21 * expon2;
  81.     Veloc12 = u11 * u12 * expon1 + u21 * u22 * expon2;
  82.     Veloc21 = u11 * u12 * expon1 + u21 * u22 * expon2;
  83.     Veloc22 = u12 * u12 * expon1 + u22 * u22 * expon2;
  84.  
  85.  /* Compute matrix M(T) for theta  */
  86.     Matrix11 = -(u11 * u11 * (expon1 - 1) / s1 + u21 * u21 * (expon2 - 1) / s2);
  87.     Matrix12 = -(u11 * u12 * (expon1 - 1) / s1 + u21 * u22 * (expon2 - 1) / s2);
  88.     Matrix21 = -(u11 * u12 * (expon1 - 1) / s1 + u21 * u22 * (expon2 - 1) / s2);
  89.     Matrix22 = -(u12 * u12 * (expon1 - 1) / s1 + u22 * u22 * (expon2 - 1) / s2);
  90.  /* miscellaneous constants -- rho is forcing and is only in iterated map */
  91.     cc1 = length1 / inertia1;
  92.     cc2 = length2 / inertia2;
  93. }
  94.  
  95.  
  96. dblrotor() {
  97.     double  moduloAB();
  98.     int     i;
  99.     int     base;
  100.     double  y_temp[20],
  101.             yval0,            /* y0, y1 are global */
  102.             yval1;
  103.     yval0 = y[0] + Matrix11 * y[2] + Matrix12 * y[3];
  104.     y_temp[0] = yval0;
  105.     yval1 = y[1] + Matrix21 * y[2] + Matrix22 * y[3];
  106.     y_temp[1] = yval1;
  107.     y_temp[2] = Veloc11 * y[2] + Veloc12 * y[3] + rho * cc1 * sin(yval0);
  108.     y_temp[3] = Veloc21 * y[2] + Veloc22 * y[3] + rho * cc2 * sin(yval1);
  109.  
  110.     for(i = 0; i < num_lyap; i++) {
  111.         base = lyapzero + vec_dim * i;
  112.         y_temp[base]
  113.             = y[base] + Matrix11 * y[base + 2] + Matrix12 * y[base + 3];
  114.         y_temp[base + 1]
  115.             = y[base + 1] + Matrix21 * y[base + 2] + Matrix22 * y[base + 3];
  116.         y_temp[base + 2]
  117.             = Veloc11 * y[base + 2] + Veloc12 * y[base + 3]
  118.             + rho * cc1 * cos(yval0) * (y[base]
  119.             + Matrix11 * y[base + 2] + Matrix12 * y[base + 3]);
  120.         y_temp[base + 3]
  121.             = Veloc21 * y[base + 2] + Veloc22 * y[base + 3]
  122.             + rho * cc2 * cos(yval1) * (y[base + 1] + 
  123.             Matrix21 * y[base + 2] + Matrix22 * y[base + 3]);
  124.     }
  125.     for(i = 0; i < lyapzero + vec_dim * num_lyap; i++)
  126.         y[i] = y_temp[i];
  127.  
  128.  /* now bring y[0] and y[1] into -pi,pi etc. */
  129.     y[0] = moduloAB(y[0], -pi, pi);
  130.     y[1] = moduloAB(y[1], -pi, pi);
  131. }
  132.  
  133. initDblRotor() {
  134.     extern int      start33 ();
  135.     extern int      dblrotor();
  136.  
  137.     vec_dim = 4;        /* the dimension of the Lyapunov vectors =
  138.                    phase space dim */
  139.     num_lyap = 0;        /* the number of Lyapunov numbers to be
  140.                    computed <= vec_dim */
  141.     lyapzero = 4;        /* y[lyapzero] is the zeroth coord of the
  142.                    zeroth lyapunov vector */
  143.     init_process = start33;
  144.     dim = lyapzero + num_lyap * vec_dim;
  145.  /* needed for rungekutta */
  146.  
  147.     X_upper = pi;        /* x scale */
  148.     X_lower = -pi;
  149.     Y_upper = pi;        /* y scale */
  150.     Y_lower = -pi;
  151.  /*  mass1,2 = C1,C2;     length2 = C4; length1=C4*sqrt(C2/ (C1+C2));
  152.     damping1,2 = C5,C6         hit time=C7   */
  153.     C1 = 1.;
  154.     C2 = 1.;
  155.     C4 = 1.;
  156.     C5 = 1.;
  157.     C6 = 1.;
  158.     C7 = 1.;
  159.  
  160.     rho = 9.;        /* forcing strength */
  161.     map = dblrotor;
  162. }
  163.